arXiv:1505.05541v2 [math.NA] 18Jun2015 


QM/MM Methods for Crystalline Defects. 
Part 1: Locality of the Tight Binding Model 

Huajie Chen and Christoph Ortner* 


Abstract 

The tight binding model is a minimal electronic structure model for molecular 
modelling and simulation. We show that the total energy in this model can be 
decomposed into site energies, that is, into contributions from each atomic site 
whose influence on their environment decays exponentially. This result lays the 
foundation for a rigorous analysis of QM/MM coupling schemes. 


1 Introduction 


QM/MM coupling methods are a class of multi-scale schemes in which a quantum mechan¬ 
ical (QM) simulation is “embedded” in a larger molecular mechanics (MM) simulation. 
Due to the high computational cost of QM models, schemes of this type have become an 
indispensable tool in many scientific disciplines [31 EH ED EU EQl Eg. The present work 
is the first part in a series establishing the mathematical foundation of QM/MM schemes 
in the context of materials modelling. 

It is pointed out in [ED that a minimal requirement for a QM model to be suitable for 
QM/MM coupling is the strong locality of forces. 


dfn 

dym 


0 


“sufficiently rapidly” as t 0, 


( 1 . 1 ) 


where = \yn — ym\ and fn denotes the force acting on an atom at position yn within a 
collection of nuclei at positions {yi} C The condition (1.1) is called strong locality to 
set it apart from the weaker condition of locality of the density matrix, which is already 
well understood (see e.g. [31 Eg and ^1.2[ ) . 

To study (1.1) we take a general tight binding method at finite electronic temperature 
as a model problem. We prove an even stronger condition than (1.1), strong energy locality. 
Given a hnite collection of nuclei y = {ye}, we decompose the total energy E = E{y) into 


( 1 - 2 ) 
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where the site energies E£{y) are local in the sense that 


dEeiy) 

< p-iren 

d^Eeiy) 

dyn 


dyndym 


< p-l{ren+ri m ) 


(1.3) 

While the specihc 
, the locality resnlt 

(|1.3l) is, to the best of onr knowledge, new. We are only aware of one analogons resnlt. 


for some 7 > 0 , and analogous results for higher order derivatives, 
form of the decomposition we employ (cf. 


2.3) is well-known 


for the Thomas-Fermi-von Weizsacker model 

This locality result has a range of consequences, such as: (1) It provides a strong 
theoretical justihcation for the concept of an interatomic potential. (2) From a purely 
analytical point of view, there is little (if any) distinction between the tight binding model 
and interatomic potential models. This means that we can apply many of the analytical 
tools developed for interatomic potential models, for example: (3) We can rigorously 
formulate and analyze models of defects in inhnite crystalline solids [27] . (4) We can 
extend the construction and analysis of atomistic/continuum multi-scale schemes. In 
particular, the Cauchy-Born continuum limit analysis [HT] can be directly applied without 
additional work. 

(5) Our main motivation, however, is to formulate and analyze new QM/MM coupling 
schemes for crystal defects. In this endeavour we build on the successful theory of atom¬ 
istic/continuum coupling |l 6 ], employing the tools and language of numerical analysis. 
The key idea is that, due to (1.2) and (|1.3[), the total energy can be approximated by 


E{y) 


5^ Et(v) + Et{y), 

teQM teMM 


where Ei is not an off-the shelf site potential (Lennard-Jones, EAM, ...) as in previous 
works on QM/MM coupling, but instead is a controlled approximation to E^. We show 
in the companion papers [HI [19] that this approach yields new QM/MM schemes (both 
energy-based and force-based) with rigorous rates of convergence in terms of the QM core 
region size. 


1.1 Outline 


In Section we focus on finite systems. We hrst present a thorough discussion of real- 
space tight binding models and then establish the results (1.2) and (1.3) in this context. 
In Section]^ we then extend the dehnition of the site energy as well as the locality results 
to inhnite systems (with an eye to crystal lattices) via a limiting procedure. In Section 
we briehy present two applications of the locality results, both in preparation for Parts 
2 and 3 of this series: We extend the crystal defect model and the convergence analysis 
for a truncation scheme from [27] to the tight binding model. Finally, in Section we 
present some preliminary numerical tests illustrating our analytical results. 


1.2 Further Remarks 

Tight binding model. Tight binding models are minimalistic quantum mechanics type 
molecular models, used to investigate and predict properties of molecules and materials 
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in condensed phases. Both in terms of accuracy and computational cost they are situated 
between accurate but computationally expensive ab initio methods and fast but limited 
empirical methods. While tight binding models are interesting in their own right, they 
also serve as a convenient toy model for more accurate electronic structure models such 
as (Kohn-Sham) density functional theory. 

The study of defects in crystals is a held to which tight binding is well suited, as it is 
frequently the case that the deviations from ideal bonding are large enough that empirical 
potentials are not sufficiently accurate, but the system size required to isolate the defect 
(e.g. for dislocations or cracks) makes the use of ab initio calculations challenging. A 
number of studies have been carried out in which tight binding is applied to simulations 
of crystal defects, see e.g. |10l ESI EHl l56] . 

Weak versus strong locality. By weak locality we mean that the electron density 
matrix has exponentially fast off-diagonal decay. In the context of tight binding, this 
means 


[r 


< p-irmr, 


(1.4) 


see (2.23) for the dehnition of the tight binding density matrix r(j/)). In physics, this 

which states that the electron 


is often described by the term “nearsightedness’ 
properties of insulators and metals at hnite temperature do not depend on perturbations 
at distant regions. This property has, e.g., been exploited to create linear scaling electronic 
structure algorithms [I11Z11331EZ]. 

However, the weak locality is not enough to validate a hybrid QM/MM approach 


We need the stronger locality condition (1.1) to guarantee that the QM region is not 


affected by the classical particles, and moreover that the forces in the QM region can be 
computed to high accuracy by only considering a small QM neighbourhood. The decay 


rate in equation (1.3) then gives a guide to how large the QM region needs to be (see also 


Theorem 4.3 and [HI IT9]). 

Thermodynamic limit. Thermodynamic limit problems (inhnite body limit), 
related to our analysis in Section have been studied at great length in the analysis 
literature. The monograph na gives an extensive account of the major contributions and 
also presents the thermodynamic limit problem for the Thomas-Fermi-von Weizsacker 
(TFW) model for perfect crystals. The thermodynamic limit of the reduced Hartree- 
Fock (rHF) is studied in (TB] for perfect crystals. This literature also contains many 
results on the modelling of local defects in crystals in the framework of the TFW and 
rHF models, see e.g. PEiEiiininiiiaiiaEi]. 

These discussions are restricted to the case where the nuclei are hxed on a periodic 
lattice (or with a given local defect). Leaving the positions of the nuclei free is also a case 
of great physical and mathematical interest. Motivated by [27] we present such a model 
in § 4.1, but postpone a complete analysis to d- 

A related problem is the continuum limit of quantum models. The TFW and rHF 
models are studied in [5l|T3] where it is shown that, in the continuum limit, the difference 
between the energies of the atomistic and continuum models obtained using the Cauchy- 
Born rule tends to zero. The tight binding and Kohn-Sham models are studied in a 
series of papers [231 [211123112S] , which establish the extension of the Cauchy-Born rule 
for smoothly deformed crystals. Our locality result yields an immediate extension of the 
analysis of the Cauchy-Born model in the molecular mechanics case |51j . 
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1.3 Notation 


The symbol {■, ■) denotes an abstract duality pairing between a Banach space and its dual. 
The symbol | ■ | normally denotes the Euclidean or Frobenius norm, while || ■ || denotes 
an operator norm. The constant C is a generic positive constant that may change from 
one line of an estimate to the next. When estimating rates of decay or convergence, C 
will always remain independent of the system size, of lattice position or of test functions. 
The dependencies of C will normally be clear from the context or stated explicitly. 


1.4 List of assumptions 

Our analysis requires a number of assumptions on the tight binding model or the under¬ 
lying atomistic geometry. For the reader’s convenience we list these with page references 
and brief summaries: 


L 

H.tb 

H.loc 

H.sym 

F 

U 

H.emb 

D 


P 

P 

P 

P 

P 

P 

P- 

P- 


15 


20 


uniform non-interpenetration 
locality of Hamiltonian 
locality of Hamiltonian derivatives 
symmetries of Hamiltonian 
conhguration independent distribution 
locality of the repulsive potential 

connection between the Hamiltonians of two embedded systems 
homogeneity of the reference conhguration outside a defect core 


2 Tight binding model for finite systems 

We begin by formulating a general tight binding model for a hnite system with N atoms. 
Let Aat be an index set with = N. An atomic conhguration is described by a map 
1 / : Aat —)■ M'’* with d G N denoting the space dimension. (We admit d ^ 3 mostly for 
the sake of mathematical generality; e.g., this allows us to formulate simplihed in-plane 
or anti-plane models.) 

We say that the map y is a proper configuration if the atoms do not accumulate: 

L. 3 m > 0 such that \y{f) — y{k)\ > nr V £, A; G Km- 

Let C denote the subset of all y G satisfying L. 


2.1 The Hamiltonian matrix 

In the tight binding formalism one constructs a Hamiltonian matrix Ti in an “atomic-like 
basis set” {^^(r - 

=[ - 2/(^))H(2/)0fc/3(r-|/(/c)) dr, (2.1) 

V / ek J^d 

where S is a small collection of the atomic orbitals per atom (with maximum size ns), and 
the exact many-body Hamiltonian operator "H is usually replaced by a parametrized one. 
The entries of the Hamiltonian matrix "H depend on the atomic species, atomic orbitals 
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and on the confignration of nnclei. In practice, they are often described by empirical 
fnnctions {empirical tight binding) which have been calibrated nsing experimental results 
or results from hrst principle calculations. 

In either case, we can write the Hamiltonian matrix elements as 

(«(»))“'’ = (a), (2,2) 

where ^ M are functions depending on i, k, a and f3. 

The orbital indices a, f3 do not bring any additional insight into the problem we are 
studying, while at the same time complicating the notation. Therefore, we ignore the 
indices a, jS, which is equivalent to assuming that there is one atomic orbital for each 
atomic site (nE:=l). The Hamiltonian matrix elements then simply become 

(n{y)) = hk{y), (2.3) 

\ / ek 


All our results can be generalized to cases with ris 
required modihcation is outlined in [Appendix A 


> 1 without difficulty. The only 


We make the following standing assumptions on the functions hik{y), which we justify 
below in Remark 2T and Examples MM Briefly, these assumptions are consistent 


with most tight binding models, with the only exception that we assume that Coulomb 
interactions are screened. 


H.tb. There exist positive constants ho and 70 such that, for any y G , 

|h^fc(l/)| < Wi,keAN. 

H.loc. There exists n > 4 such that h^k G Moreover, there exist positive 

constants hj and 7 j for 1 < j < n, such that 


(2.4) 


d^htkiy) 




- ■ ■ d[y{mj)]i^ 
with mi, • • • , mj G Atv and 1 < R,..., < d. 

H.sym. (i) (Isometry invariance) If G V({( and g : 


N 


( 2 . 6 ) 




is an isometry, then 


hek{y) = hik{g{y)) i,k e A^. 


( 2 . 6 ) 


(ii) (Permutation invariance) If ?/ G V(^ and ^iAtv— tAivisa permutation (relabelling) 
of Ajv, then 


hgkiy) — hg-i(£)g-i(fc)(|/ o Q) y i^k e An- 


(2.7) 


Remark 2.1. (i) Condition (2.4) indicates that all the matrix elements are bounded by 
ho, which is independent of the system size. This is reasonable under the assumption L 
and that the number of atomic orbitals per atom in H remains bounded as the number of 
atoms N increases. 
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(ii) The condition H.tb postulates exponential decay of the matrix elements with re¬ 
spect to the nuclei distance \y{t) — y{k)\. This is true in all tight binding models; as a 
matter of fact, most formulations employ a finite cut-off (zero matrix elements beyond a 
finite range of internuclear distance). 

(Hi) When j = 1 in H.loc, the condition (2.5) becomes 


dhkiy) 


d[y{m)]i 


< W i,ke An 


( 2 . 8 ) 


with 1 < i < d. This states that there is no long-range interactions in the models, so that 
the dependence of the Hamiltonian matrix elements hik{y) on site m decays exponentially 
fast to zero. This assumption is reasonable if one assumes that Coulomb interactions are 
screened. 

(iv) In most tight binding models, the atomic orbitals are not orthogonal, which gives 
rise to an overlap matrix 



yiWtuir - y{k))dr. 


(2.9) 


(In empirical tight binding models, AA may again be given in functional form.) 

On transforming the Hamiltonian matrix from a non-orthogonal to an orthogonal basis 
by taking the transformed Hamiltonian 

n = 


we obtain again the identity as overlap matrix. Moreover, following the arguments in m 
it is easy to see that, if AA has an exponential decay property analogous to (2.4), then so 
does AA~^I‘^. Thus, we see that the decay properties in H.tb and H.loc are not lost by 
this transformation and we can, without loss of generality, ignore the overlap matrix. 

(v) We have opted to work with an isolated system, however, it would be egually possible 
to employ periodic boundary conditions. In this case, tight binding models employ Bloch 
sums to take into account the periodic images; see, e.g., the Slater-Koster formalism 

Our entire analysis can be easily adapted to this case as well, and the resulting 
thermodynamic limit model would be identical to the one we obtain. 

(vi) H.sym (i), invariance of the Hamiltonian under isometries of deformed space, is 
true in the absence of an external (electric or magnetic) field. E.g., with g{x) = x -\- c, 
with some c G (i) implies translation invariance of the Hamiltonian. H.sym (ii) 
indicates that all atoms of the system belong to the same species so that the relabelling of 
the indices only gives rise to permutations of the rows and columns of the Hamiltonian. 

The symmetry assumptions H.sym are natural and represent no restriction of gen¬ 
erality. We reguire them to establish analogous symmetries in the site energies that we 
define in §\EE We remark, however, that H.sym must be modified for multiple atomic 


orbitals per site; see Appendix A 


Example 2.1. Many tight binding models use the ‘two-centre approximation’ fSBf, as¬ 
suming that hikiy) depends only on the vector between two atoms y{£) and y{k). If we 
only take into accounts the nearest neighbour interactions, then the Hamiltonian matrix 
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elements of such models are given by: 

( Oi if i = k; 

= < bik{y{i) — y{k)) if y{tj is the nearest neighbor of y{k)] ( 2 . 10 ) 

0 otherwise, 

where ai are constants and b^k are smooth functions. We observe that all our assumptions 
in H.tb and H.loc are trivially satisfied for this simple but common model. 

Example 2.2. The Hamiltonian of a reduced Hartree-Fock model with the Yukawa po¬ 
tential 13^ is 


«(») 


y„(.-!/(«)) + 

£eAjv 





— x)dx, 


( 2 . 11 ) 


where p is assumed to be a fixed electron density, and Ym is the Yukawa kernel with 
parameter m > 0: 



if d 

g—m|x| cosht^^ 

if d 


if d 


( 2 . 12 ) 


Note that both Ym and its derivatives decay to 0 exponentially fast. If the basis functions 
{4>ea}eeAN,aeE are localized, i.e. the atomic orbitals for the ith atom have compact support 
around y{i), or decay exponentially, then we have that the matrix elements generated by 


(2.1) satisfy the assumptions in H.tb and H.loc. 


As a conseqnence of our assumptions, the following lemma states that the sepctrum 
of the Hamiltonian is uniformly bounded with respect to the system size N. 

Lemma 2.1. For any satisfying L and H.tb, there exist constants q_ and a depending 
only on m, ho, 70 and d, such that, for all y G , 


(T{H{y)) C [a, a]. 

Proof. Using (2.4) and the Gersgorin’s theorem [SB], we have 


(2.13) 


|A,| < max I \hei{y) \ + ^ \hek{y)\ < ho 




max 

teAiv 




3-7oh(d-y(UI 


(2.14) 


for any i. This together with the assumption L implies that 

Cdho 


|Aj| < ho max IE g-7omh(d-y(UI I < 




(m7o)‘ 


(2.16) 


for any i, where Cd is a constant depending only on the dimension d. 


□ 
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2.2 Band energy 

The total energy of of a configuration y G 3^^ is written as the sum of band energy and 
repulsive energy, 


(2.16) 

which we dehne as follows. For simplicity of notation, we will write E = throughout 

this paper. 

Given a deformation ?/ G Vm, the associated Hamiltonian matrix and its eigen¬ 

values eg and eigenvectors s = 1,..., N (allowing for multiplicity), 

'H{y)il)s = s = 1,2, ■ ■ ■ ,iV, 

the band energy of the system is dehned by 

N N 

= (2.18) 

S=1 

where / depends on the physical context. For example, at hnite electronic temperature, 
/ is the is Fermi-Dirac function 

/(£) = (^1 + , (2.19) 

and yU is a hxed chemical potential (more on that below). In the zero-temperature limit, 
/ becomes a step function. In practical simulation of conductors, / is often a smearing 
function (i.e., a numerical parameter) to ensure numerical stability (see, e.g. [30l HU ITT] . 

In the present work, we shall not be too concerned about the origin of f, but simply 
accept it as a model parameter. Our analysis can be carried out whenever / is analytic 
(e.g., the Fermi-Dirac distribution) or, in insulators (systems with band gap at ju) also 
when / is a step function. For the sake of a unihed presentation we shall only present the 
hrst case, but it will be immediately apparent how to treat insulators as well. Thus we 
shall assume for the remainder of the paper that 

F. / is a conhguration independent analytic function in an open neighbourhood 
Df C C of [a, a]; cf. Lemma 

Remark 2.2. (i) The qualifier “configuration independent" in F essentially rephrases the 
assumption that the chemical potential y is independent of the configuration y. This is 
false in general, but a reasonable assumption in our context since, in the next section, we 
shall consider limits of finite bodies in the form of lattices that are only locally distorted 
by defects. It is well-known (though we are unaware of a rigorous proof) that the limiting 
potential y is indeed configuration independent, but is only a function of the far-field 
homogeneous lattice state. 

(a) We note, though, that there is a simple model in which the Fermi-level is indeed in¬ 
dependent of the configuration. Consider a single-species two-centre approximation where 
hek{y) = h{\y{i) — y{k)\) and huiy) = ao- Then it is easy to see that the spectrum is 
symmetric about Qq and hence the Fermi level is always /i = Oq. 


2.1 


(2.17) 



The repulsive component of the energy is empirical and in most of the cases is simply 
described by a pair potential interaction 


MyiO-y(k)), ( 2 . 20 ) 

where Uek is an empirical repulsive energy acting between atoms on y{i) and y{k). For 
future reference, we rewrite this in site-energy form, 

E’”'{y)^Y.E?^^y). £','•"(!/) = ^ E U„{y(e)-y(k)), (2.21) 

fc^AjY, 

and we shall assume throughout that 

U. Uik G \ B-m) and there exist cu-, yu > 0 such that 

|V-^f4fc(r)| < Ci7exp(-7(7|r|) V £, A; G Aat (2.22) 


for 0 < j < n. 

In most of our analysis we shall only be concerned with the band energy E, and have 
added mostly for the sake of completeness. The pair interaction in may be 
replaced with an arbitrary short-ranged interatomic potential. 


2.2.1 Representation via contour integrals 

Our analysis of the locality of interaction generated by the tight binding model builds 
on a representation of E in terms of contour integrals. The main issue is to represent 
the electronic density matrix as an operator-valued function of the Hamiltonian. This 
technique has been used in quantum chemistry, for example [211131] for tight binding and 
nn EH [321 ESI EZI for density functional theory. 

We begin by dehning, for any proper conhguration y G the electronic density 
matrix (or simply, density matrix) of the system. 


r( 2 /) = = fi'Hiy)). 


(2.23) 


The band energy can then equivalently be written as 

E{y) = TT{niy)T{y)) = Tr(H( 2 /)/(H(|/))) = Tr(f(H( 2 /))). 


(2.24) 


Lemma 2.1 and F imply that we can hnd a bounded contour 'if C Df, circling all the 
eigenvalues Eg on the real axis(see Figure [2T|, and satishes 


mm 


{dist('^,cr('H(?/))),dist('^,s(f))} > h, 
with a constant h > 0 that is independent of y or of N. Let 

^z{y) ■= (y-iy) -ziy 


(2.25) 
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Figure 2.1: A schematic plot of the dumbbell-shaped Cauchy contour 


denote the resolvent of then 

= 

which implies that 

E(.;ti) = 


27ii 


f{z)j^,{y) dz, 


(2.26) 


27ii 


^{z)i:i(M^{y)] dz. 


(2.27) 


It is already clear from (2.27) that the locality of the resolvents will play an important 


role in our analysis. Hence, we prove a decay estimate in the next lemma. 


Lemma 2.2. Let l-i{y),y G be a tight binding Hamiltonian of the form (2.3) and ‘to 


a contour satisfying (2.25). //L and H.tb are satisfied, then there exist constants 7 r > 0 


and Cr, independent of y or N, such that 


Ik 


(2.28) 


Proof. This proof relies on the arguments provided by |21] and a Combes-Thomas type 

estimate I2D1- 

For ko and 7r > 0, let B G 


Bik — 


From this dehnition, we have 


g7r|j/F)-2/(fco)|^ if £ = fc; 

0 , otherwise. 


(2.29) 
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Assumptions L and H.tb yield 


-nW^ < sup V - l) 


( 


< ho sup 

eeAff 


\ 




:,-{io-'yi)\yii)-y{k)\ 


fceAjv, 

\ \y(i)-yik)\>R 


+ 


E 


3-70l?/(d-?/(fc)l 1 A7r-R 


/cGAtv, 

\y{i)-y{k)\<R 


< C'(^e-5(70-7dH + g7rR_ 


1 ) 


(2.30) 


for any 7r < 70/2 and R > 0, where (7 is a constant depending only on Jiq, d, 70, 7r, 
and m. For any e > 0, we can choose R sufficiently large and then 71 . sufficiently small 
(depending on R) such that \\BRB~^ — R\\oo < e. We note that the choice of R and 
7 i. does not depend on the system size N but only on e and the constants Jiq, 70 , m. 
Similarly, we have the same bound for \\BRB~^ — R\\i. Using interpolation, we get the 
same bound for \\BRB~^ — R\\ 2 . 

Note that 


B{z-H)-^B-^ = {z-BRB-^Y^ 

= {z--{BHB-^-H){z-H)-^)-^. 


(2.31) 


Since (2.25) implies \\{z — R) ^||^(i 2 )] < 1/f), we can choose R and 7 r such that z — BRB 


-1 


is invertible and 


\\Biz-n)-^B-Y^^i2^<-. 


Using \[B{z-n)-^B-Yk\ < \\B{z-n)~^B-^\Y(^i2)<2/-0 and 

= [B{z-ny^B-^] 

and consequently, 

— I < ^g-7r(hW-y(*:o)|-|y(fc)-y(fco)l)_ 

I L J 1 ^ 

Taking ko = k, we obtain the stated exponential decay estimate. 


£k 


2 


(2.32) 

□ 


2.3 Site energy 


Since the tight binding Hamiltonian (2.1) is given in terms of an atomic-like basis set, 
we can distribute the energy among atomic sites. This is a well-known idea, which has 
been used for constructing interatomic potentials based on the bond-order concept (see 

e.g. [2Hl[2niES])- 


Noting that ||' 0 s||r 2 = 1 for all s, we have 


N 

S = 1 


N 


N 


s=i £eAjv £eAjv s=i 
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That is, we have obtained the decomposition of the band energy 

B(») = E Ei{y), where (2.33) 

Ei{y) = ^f{es)es[iJs]] = (2.34) 

s s 


and we call E^^y) the site energy. 

When the atomic orbitals are not orthogonal we slightly modify the dehnition of site 


energy slightly for compntational efficiency; see Appendix B for detailed discnssions. 
Such a decomposition is useful, for example, since classical interatomic potentials 


almost always decompose the total energy in such a way, hence the relation (2.33) can 


be used to establish a bridge between ab initio models and empirical interaction laws 
|2H]- For our own purpose, the decomposition will (1) yield a relatively straightforward 
thermodynamic limit argument to dehne and analyze variational problems on the inhnite 
lattice along the lines of 1271, and (2) provide a starting point for the construction and 
analysis of concurrent multi-scale schemes hybrid models, which we will pursue in two 
companion papers [T8l [T9] . 

Our aim, next, is to establish locality of E^. From now on, for the sake of readability, 
we drop the argument {y) in 'H{y), and 'H^mniv) whenever convenient and 

possible without confusion and in addition write r^n = \y{.fn) — y{.n)\. Let be the N 
dimensional canonical basis vector, then we obtain from (|2.17[) that 


Ee{y) = \bl’s]ef = 

s s s 


and employing (2.26) we arrive at 


Ei{y) = 


271 i 


Kz)[^z]u dz. 


(2.35) 


We can now calculate the hrst and second derivatives of E(^{y) based on (2.35) and the 
regularity assumption in H.loc: 


dE,{y) _ 1 j 


d[y{m)\i 2'Ki 


u 


dz, and 


(2.36) 


d^E,{y) 


= F- 


d[y{m)]id[y{n)]j 2'Ki 


^z -^z 

-^z 


dz. (2.37) 


ie 
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We also have higher order derivatives of the site energy for n <n\ 
d'^Eiiy) 1 


■ ■d[y{mn)]i^ 27ii 


m 


E E E (-1)' 

. 1=1 iiH- ^ji=n^hr-- Ji 

> mi, -- ,mn 




dz, (2.38) 


ee 


where Q‘f'^ : ^ 


qW[ 2 /]( 0 i, • • •, 0 .) = ^z{y) n 

i=i 


(2.39) 


is a well-defined linear map for any 2 ; satisfying (2.25) and is the mnltiset 

permntation of mi, - - - , mn- 


2.4 Properties of the site energy 


In order for to be a “true” site energy it must satisfy certain properties: locality, 
permutation invariance and isometry invariance. We establish these next. 

First, we establish the locality of the site energy and its derivatives. We remark that 
in this result it is important that we are keeping /i fixed. Admitting y-dependent /i would 
introduce a small amount of non-locality in the site-energies, but it is reasonable to expect 
that this vanishes in the thermodynamic limit. 


Lemma 2.3 (Locality). // L, H.tb, H.loc, F are satisfied then, for 1 < j < u, there 
exist positive constants Cj and rjj such that for any £ G A^r, 


d^E,{y) 

d[y{,mi)]i^ - - - d[y{mj)]i^ 




I ,ij <d. (2.40) 


Proof. We will only give the explicit proofs for j = 1,2, the cases j > 2 being analogous 
(but tedious). 

For j = 1, we have from Lemma [2. 2 [and the assumptions L, H.loc that 




ee 


^ [^z] ( [T-L,m]i [^z] 

(^iiGsAjv 


< c^hi E min{7r ,71} 


^1,-^2SAjv 


= c; 


'M E g min'[7r,7i} A 

V (J, c A / 


' ^iSA^ 


(2.41) 


where C depends only on hi, Cr, 7r, 7i, tn and d. Together with (2.36) and F, this leads to 

dEeiy) 


d[y{m)], 


< for 1 < i < d. 


(2.42) 
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For j = 2, we employ Lemma 2.2 and L, H.loc to estimate the three term arising in 


the expression (2.37) of the site Hessian, nsing similar compntations as (2.41): 


< c^hf 


^1i-^2/3/4SAjv 




It 


^z [n. 


Hj 

„ 2 ] 


ll 


< C^h2 E g- min{7r,72}(r«^+r^^„,+r^j„+r^2m+»’C2n+’'(!2'!) 

^iifesAjv 

< (^g-5 min{7r,72}(|y(£)-y(m)| + |?;(€)-?;(n)|) 


Inserting these three estimates into (2.37) together with F yields the desired resnlt, 


^ ^ (j^^-'n 2 {\yW-yim)\+\y{e)-y{n)\j for 1 < i j < d 


□ 


d[y{m)]id[y{n)]j 

The next lemma snmmarises the conseqnences of H.sym. 

Lemma 2.4 (Symmetries). Let y G . Assume that H.sym is satisfied. 

(i) (Isometry invariance) if g ■.'Mfi is an isometry, then Efiy) = Ei{g{y)); 

(a) (Permntation invariance) if Q : A ^ A is a permutation (relabelling) of A, then 
Efiy) = Eg-^ifiyoQ). 

Proof, (i) Let y' = g{y). Since g is an isometry, we have that y' G The assnmption 


H.sym (i) implies hmniy) = hmniyfi^ which together with (2.35) yields 
1 


Efiy) = 


27ri 


f{z) [^z{y)]u dz = 


2TTi 


f{z) [^z{y')]u dz = Efiy'). (2.43) 


(a) Let y' = yoQ. The assnmption H.sym (ii) implies hmniy) = hg-i(^rn)g-^(n){y'), which 
together with (2.35) leads to Ei{y) = Eg-i «)(!/')■ □ 


3 Pointwise thermodynamic limit 

Onr aim in this section is to give a meaning to energy in the inhnite body limit (“thermo¬ 
dynamic limit”). The notion of site energy makes this relatively straightforward: we will 
prove that, hxing a site i, and “growing” the material aronnd it to an inhnite body yields 
a well-dehned site energy fnnctional Ei for an inhnite body. Total energy in an inhnite 
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body is of course ill-defined, but using the site energies it then becomes straightforward 
to consider energy-differences; cf. § 

We need the following additional assumption, connecting the Hamiltonians for growing 
index-sets, in our analysis. 


H.emb. Let : Ajv — >■ y '■ A^r U {x'} — ?• be two conhgurations satisfying L, and 

hfi^{y^), hek{y) be the corresponding Hamiltonian matrix elements of these two 
conhgurations satisfying H.loc. If y^{i) = y{i) for any i G A^, then for 0 < j < n — 1, 


d[y{mi)]i^ ■ ■ ■ d[y{mj)]i^ 


lim 

| i ;( x ')|^’00 


_ d^hekjy) _ 

d[y{mi)]i, ■ ■ ■ d[y{mj)]i^ 


i,k e An 


(3.1) 


with mi, • • • , mj G A^v and 1 < A, • • • ,ij < d. 


Remark 3.1. Intuitively, H.emb states that, if one atom y{x') in the system is moved 
to infinity, then the Hamiltonian matrix elements for the remaining subsystem {y{i) \ I G 
An} do not depend on y{x') anymore. 

At first glance, this appears to be a conseguence o/ H.tb and H.loc. The reason we 
have to formulate it as a separate assumption is to make a connection between the Hamil¬ 
tonians for An and Aat U {x'}. More generally, in Lemma 3.1, we obtain an analogous 
connection between the Hamiltonians for any two systems An, Am with An C Am- 


Let A be a countable index set (or, reference conhguration), then we consider deformed 
configurations belonging to the class 

Vm(A) := {y : A^ j/|a^ C for any hnite An C A} 


= {?/ : A —>■ \y{i) — y{k)\ > m for all £, fc G A}. (3.2) 

If 1 / G Vm(A), then L is satished for any hnite subsystem A^v C A. In the following, 
whenever we assume H.tb, H.loc and H.emb for inhnite A, we mean that they are 
satished for the Hamiltonian matrices of any hnite subsystem Ajv C A. 

For a bounded domain H C M'^, we shall denote the Hamiltonian, resolvent and energy 
of the hnite subsystem contained in H, respectively, by TL^^y), and E^{y). For 

simplicity of notation, we drop the argument {y) whenever convenient. 

Theorem 3.1. Let A be countable and y G Vn,(A) be a deformation. Suppose the assump¬ 
tions F, H.tb, H.loc, H.emb and H.sym are satisfied for all finite subsystems (with 
simultaneous choice of constants), then 

(i) (existence of the thermodynamic limits) for any £ G A and for any seguence of 
convex and bounded sets H/j D BR{y{i)), R > 0 the limit 

Efiy) := hm E^^{y) (3.3) 

it^OO 

exists and is independent of the choice of sets Qr; 


(a) (regularity and locality of the limits) the limits Efiy) possess jth order partial deriva¬ 
tives with 1 < j < n — 1, and it holds that 

d^Efiy) 




< ^'-1 W)-yim,)\ fi. <d, (3.4) 


where the constants Cj and pj are same to those in Lemma 2.3: 
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(in) (isometry invariance) if g is an isometry, then Eify) = Ei{g{y)); 


(iv) (permutation invariance) if Q : A ^ A is a permutation (relabelling) of A, then 
Ei{y) = Eg-^i){yog). 


Before we prove Theorem 3.1 we establish an extension of H.emb. 


Lemma 3.1. Let Am 2 ^n, '■ Am —)■ and y^ : A^ —)■ he two configurations 

satisfying L. Assume that there exists a convex set C such that 


y^{i) = y^{i) en V £ e Atv and y^{i) G V £ G Am\Av, 


where is the complement of hi in If F, H.loc, H.ext are satisfied then, for 
0 < J < n — 1, there exist positive constants Cj, Kj, which do not depend on M, N and hi, 
such that 

\hfk{y^) - h^k{y^)\ < Co exp - Ko^dist {y^{I),Lf) + dist [y^{k),LL) (3.5) 


and 


d^hfkiy^) _ d^h^ki.y^) 

d[y{mi)]i^ ■ ■ ■ d[y{mj)]i. d[y{mi)\i^ ■ ■ ■ d[y{mj)]i. 

< Cj exp ^ — Kj ^dist [y^{i), + dist {y^{k), (3.6) 

+ ELi(I2/^W -y^i^s)] + \y^{k) -r/^(m,)|))) 


for any i, k,ms E Aj^ and 1 < ii,... ,ij < d. 


Proof. We hrst prove the case j = 0, 
a normalized vector Ui snch that y^ {i 


i.e., (3.5). For each i G Am\An, we can find 
) — ■ dist {y^{i),H) G hi. Let be a 


seqnence for each i G Am\An, snch that —)■ cxd as p —>■ cxd, then we dehne 




y^{i) if e G A^; 

y^{i) + Pe-Rf T^gAmW. 


Using H.emb with j = 0 and an elementary argnment we can indnctively choose R'f, 
snch that R^ := min£gAM\AAr t oo as /i —?• oo and 

h^,{y^) = hm hr^m- (3.7) 

11^00 


Note that the convexity of hi implies 

\y^{i) -y^(n) + ci/„| > ^(^\y^{i) -y^{n) \ +c) 
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for any n G M\N and c > 0. Using the assnmptions L and H.loc with j = 1, we have 


K(r) - = 


E 




neAM\Ajv 


dy{n) 


(1 - t)y^ + tyn ■ {y^{n) - y^{n)) dt 


< 


< 


Chi f ^ e-7i(h'“h)-y"WWi^„|+|?;"(fc)-y"(n)Wv„|) . 

n&KM\^N 


^•SAMtAjv 

where kq = 7 i /2 — £ with any £ G ( 0 , 7 i/ 2 ), and Cq is a constant depending only on £, hi, 
7 i, m and d. Note that neither kq nor cq depends on the system size M, N and f2. Note 
that 

!'>"(!/") - = lim K(v>^) - ft"®")! < 

' ' fl^OO ' ' 

which completes the proof of the case j = 0. 

With the same argnments, we can prove (3.6) for 1 < j < n — 1 by nsing the assnmp¬ 
tions H.loc with index j -|- 1 and H.emb with index J- □ 

Proof of Theorem 3.1\ (i) Withont loss of generality, we can assnme that the upper bound 
of the spectrum a < 0 (one can always shift the eigenvalues if this is not satished) and 
the contour is chosen such that it includes 0 and 


min |dist('^, a{'H{y))), dist(‘^,s(f)), dist(‘^, {0})} > c). 

Let £ G A and Bji{y{i)) C C f]', then we dehne 

\^^'^{y)\km if e VLr] 

0 otherwise for y{k),y{m) G 


(3,8) 





:= I 


km 


(3,9) 


Note that (3.8) implies that the condition (2.25) is satished for the Hamiltonian with 
the contour Moreover, the resolvent = {'hP^{y) — zl)~^ is well-dehned for any 

z G ^ and satishes the estimate in (2.28). Under the assumption F, we can observe that 
the band energy and site energies of the Hamiltonian PL^^{y) and PL^’^{y) same. 
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We have from (|2.35|), (|3.9|), L, H.tb, Lemma |3.1| and Lemma |2.2| that 


f(^) 


i£ 


<c [•*?']«. [«"' - 


< 


^ E g- 7 rr«j g-Ko (dist(y(ri),0'\0fl)+dist(y(£2),f1'Wfl)) g- 7 rr «2 


y{il),y{e2)&^' 


< C g-7rr’££i-Kodist(j;(£i),0'\nfl) j ^ q mm{'yr, ko}R 

yyiii)£n' 


(3.10) 


where the last constant C depends only on hg, Cr, Cg, 7r, nr and d, bnt is independent 
of y or R. Since (3.10) holds for any f2' 2 it follows that is a 


Canchy seqnence. The nniqneness of the limit is also an immediate conseqnence of the 
fact that fl' was arbitrary. This completes the proof of (i). 

(ii) Case j = 1: 

For i,m E A, we take R > 2\y{m) — y{i)\, and then adopt the notation in the proof 


of (i). With the expression of (2.36), we obtain by a direct calcnlation that 


dE^\y) dE^-{y) ^ ^ 




d[y{m)\i d[y{m)\i 2m 


- n^') ( [n%]^ - 


dz. (3.11) 


u 


Using L, H.tb, Lemma |3.1| and Lemma |2.2[ we can obtain from a similar argnment as in 


(3.10) that 


dE^'iy) dE^^iy) 


d[y{m)]i d[y{m)]i 


< 


(Jq~ mi°{7r,K0 All (^J+Lh)-?/(m)l) /2 


(3.12) 


where the constant C depends only on hg,Cr, cg, ci, 7r, Kg, ki, m and d. Note that 
the estimate in (3.12) can be bonnded by which does not depend on 


y. Therefore, we have that {dE^{y)/d[y{m)]i^ converge nniformly to some limit, 
which together with (i) implies that E^{y) is differentiable with respect to ?/(m) and the 
derivative is given by 


dEejy) 

d[y{m)]i R^lo d[y{m)]i 


for 1 < i < d. 


(3.13) 


Case j > 1: 
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For the second order derivatives, we can obtain from the expression (2.37) and a 
tedious calculation that 

d^Efiy) 


_ d^Ef^jy) 

d[y{m)]id[y{n)]j d[y{m)]id[y{n)]. 


= ^ 








+ 


•«?' - [^%]) «?" [«?],■«?' 






+ 


dz. 




For readability, we have omitted eight other terms in the square bracket, which have 
the same structure as the listed four terms. By estimating each term using the same 


arguments in (3.10), we can obtain 

d^Ef\y) d^Ef^iy) 


< 


d[y{m)]id[y{n)]j d[y{m)]id[y{n)]j 
which leads to the existence of d'^Ei^y)/d[y(jn)]id[y{n)]j and 
d^E,{y) d^E^^iy) 


min{7r,Ko,Ki,K2} {r+ti m +r-,„)/4^ (3.14) 


= lini 


d[y{m)\id[y{n)\j R^oo d[y{m)]id[y{n)]j 
For 2 < j < n — 1, we can obtain by similar arguments that 
d^Eiiy) ^ d^Ef^iy) 


for 1 < f, j < d. 


(3.15) 


1 < • • • , < d (3.16) 


d[y{mi)]i, ■ ■ ■ d[y{mj)]i. R^oc d[y{mi)]i, ■ ■ ■ d[y{mj)]i 

for any sequence Qr satisfying the conditions of (i). 

The locality of E( in (ii) is now an immediate consequence of Lemma 2.3 and (3.13), 
( l3T5| ), ( |3T6l ). 

(Hi) Let y' = g{y). Since g is an isometry, we have that g[BR{y{i))^ = Bn^y'^i)) for any 
d? > 0. We can obtain from H.sym (i) and Lemma 2.4[i) that 




Taking the limit d? —)• cx) of (3.17) and (i) yield E£{y) = Ek{y'). 

(iv) Similar to the proof of (Hi), (iv) is a consequence of Lemma |2^ii) and (i). 


(3.17) 


□ 


Theorem 3T states the existence of the thermodynamic limits of the site energies, as 
well as the regularity, locality and isometry/permutation invariance of the limits. In the 
following we shall always denote this limiting site energy by 77^. 
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Remark 3.2. We have only considered the band energy of the system so far. The repulsive 
energy can be incorporated into our analysis without difficulty. 


Using the expression (2.21) and the assumption U, it is easy to justify the thermody¬ 


namic limit of the repulsive site energy as well as its regularity and locality as those 
in Theorem 3.1. Moreover, the symmetry results in Theorem \3.1\ are also clearly satisfied 
with the expression (2.21) for E^^. Therefore, all we have to do is to take the total site 
energy 

= Ei + p 

and then use the existing results for Ei. For convenience and readability, we still work 
with the site band energy E^ and continue to ignore the repulsive component. 


4 Applications 

4.1 Tight-binding model for point defects 

As alluded to in the Introduction, our primary aim in understanding the locality of the TB 
model is the construction and rigorous analysis of QM/MM hybrid schemes for crystalline 
defects, along the lines of izg. The next step towards this end is a rigorous dehnition of 


a variational problem that is to be solved. Since we have shown in Theorem |3.1| that the 
total TB energy can be split into exponentially localised site energies, this is a relatively 
straightforward generalisation of the analysis in [27], which considers MM site energies 
with bounded interaction radius. 

Here, we only summarize the results, with an eye to the application we present in 
§0 For simplicity we restrict ourselves to point defects only. Complete proofs and 
generalisations to general dislocation structures are given in HZI. 

We call an index set A a point defect reference configuration if 

D. 3 Rdef > 0, A G SL((i) such that and A n B^^^^ is finite. 

While, in previous sections, we have worked with deformations y where y{i) denotes 
the deformed position of an atom indexed by i, it is now more convenient to work with 
displacements n : A —)■ u{i) = y{i) — 1. For displacements u we dehne the energy- 

difference functional 

£{u) ■.= '^[Ei{x + u) - E^{x)\, (4.1) 

£eA 

where x denotes the identity map x : A — )■ = £. Due to the exponential locali¬ 

sation of El this series converges absolutely if u has compact support, i.e., for u G 
with 

:= {u : A —)■ Mf, 3R > 0 s.t. u = const in A \ Hr}; 


cf. Theorem 4.1 


1 . 


Next, still following [27], we extend the dehnition of T to a natural energy space. For 
i & A, p E A — i := {m — i, m G A\{£}} we dehne Dpu{i) := u{i-\-p) — u{i), and moreover. 


Du{i) := {Dpu{i)) 


pSA-G 
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We think of Du{i) G (M^)^ ^ as an (infinite) finite-difference stencil. For any snch stencil 
Du{i) and 7 > 0 we define the norm 

peA-e 

which gives rise to an associated semi-norm on displacements, 

/ \ 1/2 

^ ieA ' 

All (semi-)norms || • ||£ 2,7 > 0, are eqnivalent. With these dehnitions we can now dehne 
the fnnction space, which encodes the far-held bonndary condition for displacements, 

:= {m : A ^ \\Du\\p^ < 00 }. 

We remark that is dense in W d]. 

In addition to the decay at inhnity imposed by the condition u G we also reqnire 
a variant of L, stating that atoms do not collide. Thns, onr set of admissible displacements 
becomes 



Adnim := {u G -f u{t) — m — u{m)\ > m\£ — m| V £, m G A}, 


where m is an arbitrary positive nnmber. Again, we can observe that l^'^nAdmn, is dense 
in Admm na. We remark also that, dne to the decay imposed by the condition u G 
if M G Admo then u G Adm^ for some m > 0. 

We can now state the main resnlt concerning the energy-difference fnnctional £. The 
proof is an extension of m Lemma 2.1] and will be detailed in [1^. The main new 


ingredient in this extension, as well as in Theorem 4.2 below, is qnantifying how rapidly 
the site energies approach those of a homogeneons crystal (withont defect). 

Theorem 4.1. Suppose that D is satisfied, as well as F, H.tb, H.loc, H.emb and 
H.sym for all finite subsystems with simultaneous choice of constants. 

(i) £ ■ n Admo —)■ M is well-defined by (4.1), in the sense that the series converges 
absolutely. 

(a) £ : n Admo —)■ M is continuous with respect to the || • \\i 2 semi-norm; hence, 

there exists a unique continuous extension to Admo, which we still denote by £. 

(in) £ G C"“^(Admo) in the sense of Frechet. 


In view of Theorem 4.1 the following variational problem is well-dehned: 


u G argmin |T(m), m G Admo}, 


(4.2) 


where “argmin” is nnderstood in the sense of local minimality. We are not concerned 
with existence or nniqneness of minimizers bnt only their strnctnre. This is discnssed in 
the next result, which is an extension of m Thm. 2.3] (see [T7] for the complete proof). 
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Theorem 4.2. Suppose that D is satisfied, as well as F, H.tb, H.loc, H.emb and 
H.sym for all finite subsystems with simultaneous choice of constants. If u & Admo is a 


strongly stable solution to (4.2), that is, 


3 c > 0 s.t. (^6‘^£{u)v, o') > c||i3u||^2 Vu G ) 
then there exists a constant C > 0 such that u satisfies the decay 

\Du{I%<C(l + \I\)-^ We A. 


(4.3) 


Remark 4.1. (i) The condition ( |4.3 ) is stronger than actually reguired. Indeed, it suffices 
that u is a critical point and that strong stability is satisfied only in the far-field; cf. 

Eg. (2.7)]. 

(a) Higher-order decay estimates can be proven for higher-order gradients. For ex¬ 
ample, when n > 5, then \Dp.^Dp.^u[I)\ < for \I\ sufficiently large; see [TT}/ 


for more details. These estimates will be useful in our companion papers uacw for the 
construction of highly accurate MM potentials, but are not reguired in the present work. 


4.2 Convergence of a numerical scheme 

As a reference scheme to compare our QM/MM schemes against, and also as an elementary 
demonstration of the usefulness of the locality results and of the framework of § 4A we 
present an approximation error analysis for a basic truncation scheme. 

To construct the scheme we hrst prescribe a radius R > 0 and restrict the set of 
admissible displacements to 

Admo(i?) := {u G Admo, m = 0 in A \ Br]. 

The pure Galerkin scheme ur G argmin{£^(M ),m G Admo(i?)} is analyzed in and the 
convergence rate \\Du — DuR\\i2 < jg proven. 

In our case, the energy-difference £{u) is not computable for u G Admo(i?) due to 
the inhnite interaction radius of the TB model. However, we can exploit the exponential 
localisation to truncate it. To that end, we let Rf"^^ be a buffer region width (cf. Theo¬ 
rem 4.3 and Remark 4.2), Ar := AflR^+^buf and for any : A —)■ dehne : Ar — )■ 

satisfying = v on Ar. Then, for u G Admo, we dehne the truncated energy-difference 
functional 

£r{u) := E^^{[x + uf) -E^^{x^). 

Clearly, Sr is well-dehned and Sr G C"“^(Admo) in the sense of Frechet. To formulate 
the computational scheme, Sr need only be dehned for u G Admo(i?), but for the analysis 
it will be convenient to dehne it for all u G Admo. 

The computational scheme is now given by 


Ur G argmin 


{£r{u),u G Admo(R)}. 


(4.4) 


Theorem 4.3. Suppose that D is satisfied, as well as F, H.tb, H.loc, H.emb and 
H.sym for all finite subsystems with simultaneous choice of constants. 
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Ifu is a strongly stable solution to (4.2), then there are constants C, Rq, Cbuf such that, 
for R > Rq and > Cbuflog(i?), there exists a strongly stable solution ur to (4.4) 
satisfying 


\\Du-Dur\L <cr-'^/^ 

II 11^^ 

\S{u) - Sr{ur)\ < CR 


and 


-d 


(4.5) 

(4.6) 


Proof. We closely follow the classical strategy of the analysis of hnite element methods, 
which is detailed for a setting very close to ours in [ 2 ^ in various approximation proofs. 

1. Quasi-best approximation: Following m Lemma 7.3], we can construct Tru G 
Admo(i?) such that, for any 7 > 0 and for R sufficiently large, 

\\DTru - Du\\,. < < CR-^/^ 


where Theorem |4.2| is used for the last inequality. We now £x some r > 0 such that 
Br{u) C Admrn for some m > 0. Then, for R sufficiently large, we have that Tru G 5^/2(h) 
and hence Br/ 2 {TRu) C Adm^. 

Since £ G C^(Admo(i?)), 5£ and 5‘^£ are Lipschitz continuous in Br{u) fl Admo(i?) 
with Lipschitz constants Li and L 2 , that is. 


\\6£{u) - 6£{Tru)\\ < L,\\Du - DTr{u)\U 2 < CR-'^/\ 
\\5^£{u) - 5^£{Tru)\\ < L 2 \\Du - DTr{u )\\,2 < CR-'^/\ 


and 


(4.7) 

(4.8) 


2. Stability: Using (3.14) and the facts that n = 0 and Tru = 0 outside Br, we have 
that there exists a constant 7 s, such that 

\{{6^£r{Tru) - 6^£{Tru))v,v)\ < (4.9) 

The proof of this identity is relatively straightforward but does require some details, which 
we present following the completion of the proof of the theorem. Together with (4.3) and 
(4.8) this leads to 

S^£r{Tru)v,v) 

= {S^£{u)v,v) + {{S‘^£{Tru) - S‘^£{u))v,v) + {{S‘^£r{Tru) - S‘^£{Tru))v,v) 

> {c-C{R-^/^ + e-^^^'^'''R'^))\\Dv\\%^ > ^\\Dv\\%^ Vn G Admo(i?), (4.10) 


for sufficiently large R and sufficiently large Cbuf- 

6 

that 


3. Consistency: Similarly to (|4.9|), we can derive that there exists a constant 7 c such 

(4.11) 


|(4£H(T„fi) - SS(Tru),v)\< Ce-''-6 " td-^l^\\Dv\\e. 

We also present the detailed proof of (4.11) after the proof of this theorem. 

In order to ensure that the truncation of the electronic structure (the variational crime 
committed upon replacing £ with £r) we must choose such that < 

CR-'^R^ or equivalently, . On taking logarithms, we observe that 

this is true provided that > Cbuf log R for Cbuf sufficiently large. 
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Next, employing (4.7) we obtain that 

6Sr{Tru),v) = {68r{Tru) - 6S{Tru),v) + {6S{Tru) - 6S{u),v) 

< \\Dv\\e 2 < CR-'^^^WDvW^ V n e Admo(i?) (4.12) 

for sufficiently large R and appropriate Cbuf- 

4- Application of inverse function theorem: With the stability (4.10) and consistency 


(4.12), the inverse function theorem jlHl Lemma B.l] implies the existence of ur and the 
estimate ( |4.5 ). 

5. Error in the energy: For the estimate of energy difference functional, we have 
from £ G C'^(Admo) that 

\£{ur) - £{u)\= f {6£{{1 - s)u + sur),ur- u) ds 
Jo 

j {^5£{fl — s)u + sur) — 5£{u),ur — u) ds < — Zl-u ||^2 < (4.13) 


Using (3.10) and = 0 in A\Br, there exists some constant 7 e such that 

\£{ur) — £r{ur)\ 

{Ee{x + u^)-E^^{x + u^) + Ef^{x)-Ei{x)) 

£eAni?^^^buf /2 

+ {Ee{x + u^) - Ee{x)) + ^ {Ef^{x + u^)-Ef^{x)) 


< c 


< CR^-^e 


,^£eAni?^_i_^buf /2 
d-l -7eR’=“V2 




E 


,-7e(|d-R) 


reA\B^_i_^buf /2 


(4.14) 


Using (4.13), (4.14), and possibly choosing a larger constant Cbuf, we obtain (4.6) for 
sufficiently large R. □ 


Proof of (Q. Let u := Tru and Vp := Vp{Du{^)) := Ef{x + u)-Ef{x) and W := U/. 
Using M = 0 in A \ Br we have 

{{5^£r{u) - 5‘^£{u))v,v) = Y. ((<5V/«-5V,)Z1u(£),Z1u(£)> 

^eAnSfl 

+ E ((<5V/«-5V,)Z1u(£),Z1u(£)>- {{5^V^)Dv{^),Dv{^)) 

IgArXEr ieA\AR 

=: T 1 + T 2 + T 3 . 

It is obvious from ( 3.14[ ) that 

ieAnBR peA-£ cr&A-i 


< Ce-^^'^'^'WDvW 
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with 7 = min{7r, kq, ki, K2}/4. Using n = 0 in A \ Br, we have from (3.14) that 
^ 


To < 


IGArXBr peA-£ (t£A-£ 

£+pGBji £+a£Bji 


< 


< 


£&Aii\Bii \meAnBji n&AnBji ) 

\mgAnS_R nsAnBij / 


and from (3.4) that 

T3<C J2 Y1 Y1 e-^"^^p\+\"^^\Dpv{^)\\Dc 

£&A\Ar P 6 A-^ aGA-e 

^ l+pGBji t+a£Bji 


-V 


<c Y ( H E e-’'»'l'-’"l+l'-"lA \Dv(i)\l„ 

£eA\Afl VmeAnSfl nsAnBfl / 

/ \ 1/2 


-2??2A 


)buf 


E E » 

\m£AnBji nsAnfJfl 




□ 

= 0 


772/i^ ' 

The estimates for Ti, T2 and T3 together yields (4.9) with js = niin{7,72}- 

Proof of (4.11). We continue to adopt the notation from the proof of (4.9). Using u, v 
in A \ Br we have 

{6Sr{u)-6S{u),v)= {SV,^^-SVe,Dv{e))+ - SVe, Dv{e)) 

£eAnBn 

- 5 ^ {SV,,Dv(£)) 

e&A\Aii 

=: T1 + T2 + T3. 

From (3.12) it follows that, for £ G A fl Br, — U,p| E with 

7 = min{7r, Kq, Ki}/ 2, and hence, 

TiEC e-E«+«"'-l'l+|p|)|7)^^;(£)| 

££Aji\Bji p&A—£ 

££AnBji 


sc( E^' 

' AnSfl 


1/2 
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Similarly, using u = 0 in A \ we have from (3.12) that 


T2<C 


P 6 A-^ 


< Ce-2 




E E- 


,7(K|-|^-fc|) 

keBji 

1/2 

^ ^ e’l‘l ] ||£>«|| 

i £sAfl\B/j kGBji 


1/2 / 

E E 


\ 


1/2 


e-^^P^\DpV 


11 ^^^{Kr\Br) 


< Drf-itDbuf\V2 


Finally, from 


we obtain 


T3<C E E 


e-'^^\p\\DpV 


ieAXAu peA-e 

« l+pSBn 


E E- 

i£&A\Aii k£Bji 


1/2 


/ 


3-»?i|r-fc| 


\ 


1/2 


E E 

yeeA\AR 


e-P^\p\\DpV 


peA-e 

e+psBn 


1/2 


< c 


E E ^ 


^^-ri^{\e\-R) 






The estimates for Ti, T 2 and T 3 together yields (4.11) with 7 c = min{ 7 , 7 i}/ 2 . 


□ 


Remark 4.2. The choice of buffer width R°™ is the most interesting aspect of Theo 


rem 


4-3. As expected from the exponential localisation results, we obtain that should 


be proportional to \og{R). The fact that the constant of proportionality is important 
makes an implementation difficult. At least according to our proof, if we were to choose 
Rbni _ (;;log(R) with a small c, then we would obtain a reduced convergence rate. 


Our numerical results in § |5.3| show no such dependence, which may indicate that 
our proof is in fact sub-optimal, however it is egually possible that this effect can only be 
observed for much larger system sizes than we are able to simulate. 


5 Numerical results 

We present numerical experiments to illustrate the results of the paper: (1) the locality 
of the site energies and (2) the convergence of the truncation scheme described in § 4.2 


Given the present paper is primarily concerned with the analytical foundations, we will 
show only a limited set of results, employing a highly simplihed toy model. We will present 
more comprehensive numerical results in the companion papers m [T9] . All numerical 
experiments were carried out in Julia [1]. 
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5.1 Toy Model 

The Hamiltonian matrix is given by 

Him{y) = h{\yi-ym\) where h{r) 

and /cut(r) 

with model parameters a = 2.0, ro = 1.0,rcut = 2.8. The pair potential term is set to be 
zero. 

Numerical tests suggest that a triangular lattice with 


_+ ,r<rcut, 

0, r > Tcut, 


— S 


1 1/2 \ 
0 ^/ 2 ) 


and s a scaling factor (close to 1.0) is a stable equilibrium in the sense of § 4.1 


The 2D setting and the single orbital per site make this a convenient setting for 
preliminary numerical tests. 


5.2 Locality of the site energy 


We construct two test conhgurations: (1) We “carve” a hnite lattice domain Kr = 
from the triangular lattice, and perturb each position y G Kr by a vector with 
entries equidistributed in [0, 0.1] to obtain y. (2) We obtain a second test conhguration by 
removing some random lattice sites from Kr (vacancies) and perturb the remaining posi¬ 
tions as in (1) to obtain y. We then compute the hrst and second site energy derivatives 
Eo^miy) and i?o,mn(|/) and plot them against, respectively, rom and ro™ ron- 


In the test shown in Figure [572| we chose s = 1.0 and R = 10.0, and the sites removed 
in (2) are Ari(l, 0), v4tri(0, —3), Htri(—2, 2) and Htri(2, 5). We clearly observe the predicted 
exponential decay. 


5.3 Convergence rate 


In our second numerical experiment we conhrm the prediction of Theorem 4^ We adopt 
again the model from 


tion. 


5.1 As reference conhguration we choose a di-vacancy conhgura- 
A = AtHZ2\{(0,0),(l,0)}. 


Then, for increasing radii R with associated buffer radii R 


buf 



R 

3 

4 

6 

8 

11 

Set 1 


2.1 

2.4 

2.8 

3.0 

3.4 

Set 2 

Rhnl 

1.0 

1.7 

1.7 

2.0 

2.0 

Set 3 

Rhul 

1.0 

1.0 

1.7 

1.7 

2.0 


(5.15) 


we solve the problem (4.4). In Set 1 we have chosen = 1 -|- log(i?), while in Sets 2, 3 
we have chosen smaller buffer radii to investigate the effect of these choices on the error 


in the numerical solution. 


27 












Distance [r^^+r^J 


Figure 5.2: Locality of the site energy for the tight-binding toy model described in 
Red dots denote conhguration (1), while blue crosses denote conhguration (2). 


5.1 


The computed solutions ur are compared against a high-accurac y so lution with R = 
20, = 11, which yields the convergence graphs displayed in Figure fully conhrming 

the analytical prediction. To measure errors, instead of ||Z1 ■ 11^2, we employ the equivalent 
norm ■ \\p with = {Dpu)p^±A,,;ea=i, 2 - 

We do not observe a pronounced buffer size effect. This could have a number of 
reasons, such as the fact that we are not far enough in the asymptotic regime or simply 
that the model we are employing is “too local” to observe this. We will present more 
extensive numerical results with a wider variety of TB models in [181 HH] ■ 


6 Concluding remarks 


The main purpose of this paper was to set the scene for a rigorous numerical analysis 
approach to QM/MM coupling. We have achieved this by developing a new class of 
locality results for the tight binding model. Precisely, we have shown that the total band 
energy can be decomposed into contributions from individual sites in a meaningful, i.e. 
local, way. 

This strong locality result is the basis for extending the theory of crystalline defects 
of [27], which we have hinted at in § 4.1, and carried out in detail in HU. In further 


forthcoming articles [T8l[T9] we are employing it to develop new QM/MM coupling schemes 
for crystalline defects as well as their rigorous analysis. 

A key question that remains to be investigated in the future, is whether our locality 
results extend to more accurate electronic structure models such as Kohn-Sham density 












Figure 5.3: Convergence of (4.4) with increasing domain size R, as described in § 5.3 Set 
1: full lines; Set 2: dashed lines; Set 3: dotted lines. 


functional theory. Understanding this extension is critical to take the theory we are 
developing in the present article and in [181 EH] towards materials science applications. 
However, there are many technical issues arising from the nonlinearity, the continuous 
nature, and in particular the long-range Coulomb interaction. 


Appendix A Multiple orbitals per atom 


We have assumed in § |2.1[ and throughout this paper, that there is only one atomic 
orbital for each atom (ns=l). In this case, the symmetry assumption H.sym (i) is 
natural. However, in practical calculations, there are multiple atomic-like orbitals (pia 
associated with each atomic site i. Here, a denotes both the orbital and angular quantum 
numbers of the atomic state. Many TB models employ one s-orbital |s), three p-orbitals 
{|a:), \y), |z)} and hve d-orbitals {|x|/), \yz), \zx), \x^—y‘^), |3z^ —r^)} per atom PU1I34] . 
Except for the s-orbital, all other orbitals have a spacial orientation, which means that 
the Hamiltonian matrix elements are not invariant under rotation/reflection. Therefore 
H.sym (i) is invalid and must be reformulated. 

We assume in the following that ns > 1 and {0to}i<a<n3 is the set of atomic-like 


orbitals for the site The Hamiltonian can be expressed by (2.1), 


Wfe)j 


\ 


Ik 


(r - y{^))R{y)(t)kp{Y - y{k)) dr. 


(A.l) 
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Applying an isometry g to y, we obtain 

('^{aiy))) ^ = [ 0 to(r-^(|/(^)))'H(^(?/)) 0 fc/ 3 (r-^(l/(fc))) dr. (A.2) 

V / ik J^d 

For simplicity of notation, we define 

= (j)ia(^g~\r) - y{i)^ and (r) = (/>£„(^r - 51 (?/(£))). 

We assnme that the two sets of atomic orbitals {'il^ia}i<a<ns and {v3to}i<o<ns span the 
same snbspace. This is trne for almost all tight-binding models since the set of the 
atomic orbitals always inclnde all three p—orbitals (if the p—orbital is involved) and all 
hve d—orbitlas (if the d—orbital is involved). Then, there exists an orthogonal matrix 
Qi g j^^sxns Qayi’ey- We have from (A. 2 ) that 

. f -- 

/ ‘Pia{r)'H(g{y))ipkg(r) dr 


W(9(!/)) 


ek 


Y1 Qic^'Q^' / '4’ea'{r)n{g{y))'ipky'{r) dr. (A.3) 

l<a',l3'<ns 


Since g is an isometry, it is natnral to assnme that 

'^iy)) ^ = [ ^to(r)'H(p(p))V>fc/ 3 (r) dr. 

/ £k Jmd 


(A.4) 


Let [niy) 


Ik 


W(!/) 


q/3 


£k 


denote the local Hamiltonian, then 


l<Q:,/3<n5 


we have from ( |A.3 ) and (A.4) that 

('H{g(y)]) =Q‘-(n(y)) ■ (Q*)’', (A. 6 ) 

which yields 

'^{giy)) = Q ■'^{y) ■ with Q = diag{Q\--- ,Q^} . (A. 6 ) 

Note that are orthogonal matrices, hence Q is orthogonal as well. Therefore, the 
spectrnm of 'H(p) and T-L^g^y)) are eqnivalent: 

(A.7) 


6 = 


for 1 < s < iV ■ ? 7 ,= . 


Let 




/ ^« \ / \ 

with = 


/ 


V 




£n~ 


for 1 < s < ■ n= 


/ 


be the eigenfnnction of 'H(p) corresponding to the eigenvalne e^. Then the corresponding 
eigenfnnction of T-L[g{y)) is 


= 


/ \ 


V / 


(A. 8 ) 
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Next we note that, with multiple orbitals per atom, the expression (2.34) should be 
rewritten as 


Ei{y) = E E = E E ■ 


(A.9) 


Taking into account (A.7), (A.8) and (A.9) we obtain invariance of the site energy under 
isometries. 


E,{y) = Ec(g(y)), 


(A. 10) 


To summarize, in the case of multiple orbitals, the assumption H.sym (i) should 
become 


H.sym’ (i). If ?/ G V,„ and g : ^ 

orthogonal matrices G 

(This is equivalent to H.sym (i) when = 1. 


^ is an isometry on M'^, then there exist 
nsxns for £ = 1, • • • ,77 such that (A. 61) is satished. 


Remark A.l. Slater and Roster worked out expressions such as (A.5) and (A.6) for all 
integrals between s, p and d orbitals and presented them in Table 1 of their paper 
This has been invaluable for practical calculations, see, e.g. 


Remark A.2. ITe stress again that all our assumptions and analysis in the present paper 
can be extended to the multi-orbital case without any difficulty, by taking the Hamiltonian 
as a block matrix with 


hik{y) = i'Hiy) 


ek 


(All) 


and \hik{y)\ the Frobenius norm of the submatrix. 


Appendix B Site energy with non-orthogonal orbitals 


We consider the tight binding model with non-orthogonal atomic orbitals in this appendix. 
It has been shown in Remark 2.1 (iv) that the transformed Hamiltonian is 




when the overlap matrix is not an identity matrix. Then the transformed eigenvectors of 
TLfs = become ips = and following (|2.33), the site energy is given by 




(B.l) 


Since the square root of a matrix in (B.l) introduces additional computational cost for 
the site energy computations, we modify the definition of site energy in practice by 

E,{y) = Y,K^s)[M^PsWsV (B.2) 


The following result states that the modified site energy (B.2) preserves the locality prop¬ 
erty. 
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Lemma B.l. Assume that L, H.tb, H.loc, F are satisfied, and moreover, the overlap 
matrix M. satisfy the same conditions as those in H.tb, H.loc. 

Then, for I < j < n, there exist positive constants Cj and fjj such that for any £ G A^r, 


d^Efiy) 

■ ■ ■ d[y{mj)]i^ 




1 < A, • ■ ■ fij < d. 


(B.3) 


Proof. Let S’ = f(fH) (with Sjk = XlsThe assumptions on PL and Ai 
imply that the transformed Hamiltonian PL also satishes the conditions in H.tb and 
H.loc. Using Lemma 2S and similar arguments as those in the proof of Lemma 2^ we 
have 


< C'e-^l^(T-2/(UI and 


dl 


pk 


< (jQ--f{\yin)-yU)\+\yin)-y{k)\) 


d[y{n)]i 

with some constants C and 7. Similarly, the assumptions on Ai also imply 


(B.4) 


<Ce-'^\y^^Eyik)\ and 


dM 


±1/2 

jk 


d[y{n)]i 


< (jQ-l{\y{‘n)-y{j)\ + \y{n)-y{k)\) ^ 


(B.5) 


We have from ( |B.1[ ) and ( |B.2[ ) that 

= f(c.) Y.AA'AUi’.V = 

s jk 

Ee = '^f{£s)'^Mej[i)s\j[4^s]e = 

s j 


Therefore, 

dEe 

d[y{n)]i 


E 

jk 


dM 


1/2 

u 


d[y{n)]i 




- 1/2 


(1/2 


dl 


pk 


M 


-1/2 


ke 


+ m; 


1 / 2 ; 


dM 


-1/2' 


ij -jk 


kl 


d[y{n)\i 


which together with (B.4), (B.5), and a similar argument as that in (2.41) completes the 
proof of (B.3) for j = 1. 


The proofs for 2 < j < n are similar. 


□ 


References 

[ 1 ] R. Baer and M. Head-Gordon, Sparsity of the density mtrix in kohn-sham 
density functional theory and an assessment of linear system-size scaling methods, 
Phys. Rev. Lett., 79 (1997), pp. 3962-3965. 

[2] M. Benzi, P. Boito, and N. Razouk, Decay properties of spectral projectors with 
applications to electronic structure, SIAM Review, 55 (2013), pp. 3-64. 

[3] N. Bernstein, J. Kermode, and G. Gsanyi, Hybrid atomistic simulation meth¬ 
ods for materials systems. Rep. Prog. Phys., 72 (2009), pp. 26051 1-25. 


32 























[4] J. Bezanson, a. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh 
approach to numerical computing, arXiv e-Prints, 1411.1607 (2014). 

[5] X. Blanc, C. L. Bris, and P.-L. Lions, From molecular models to continuum 
mechanics, Arch. Rat. Mech. Anal., 164 (2002), pp. 341-381. 

[6] - , On the energy of some microscopic stochastic lattices, Part I, Arch. Rat. Mech. 

Anal., 184 (2007), pp. 303-340. 

[7] D. Bowler and T. Miyazaki, 0(N) methods in electronic structure calculations. 
Rep. Progr. Phys., 75 (2012), pp. 036503-036546. 

[8] E. Ganges and C. L. Bris, Mathematical modeling of point defects in materials 
science. Math. Models Methods Appl. Sci., 23 (2013), pp. 1795-1859. 

[9] E. Ganges, A. Deleurence, and M. Lewin, A new approach to the modelling 
of local defects in crystals: The reduced Hartree-Fock case, Gommun. Math. Phys., 
281 (2008), pp. 129-177. 

[10] -, Non-perturbative embedding of local defects in crystalline materials, J. Phys.: 

Gondens. Mat., 20 (2008), pp. 294213 1-6. 

[11] E. Gances and V. Ehrlacher, Local defects are always neutral in the Thomas- 
Fermi-von Weiszdcker theory of crystals. Arch. Ration. Mech. Anal., 202 (2011), 
pp. 933-973. 

[12] E. Gances, S. Lahbabi, and M. Lewin, Mean-field models for disordered crystals, 
J. Math. Pures Apph, 100 (2013), pp. 241-274. 

[13] E. Gances and M. Lewin, The dielectric permittivity of crystals in the reduced 
Hartree-Fock approximation. Arch. Ration. Mech. Anal., 197 (2010), pp. 139-177. 

[14] E. Gances and N. Mourad, A mathematical perspective on density functional 
perturbation theory. arXiv:1405.1348. 

[15] 1. Gatto, G. L. Bris, and P.-L. Lions, The Mathematical Theory of Thermo¬ 
dynamic Limits: Thomas-Fermi Type Models, Oxford Mathematical Monographs, 
Hardcover, 1998. 

[16] -, On the thermodynamic limit for Hartree-Fock type models, Ann. 1. H. Poincare, 

An., 18 (2001), pp. 687-760. 

[17] H. Ghen, Q. Nazar, and G. Ortner, Variational problems for crystalline defects. 
in preparation. 

[18] H. Ghen and G. Ortner, Construction and analysis of energy-based QM/MM 
hybrid methods, in preparation. 

[19] -, Construction and analysis of force-based QM/MM hybrid methods, in prepara¬ 

tion. 


33 



[20] J. Combes and L. Thomas, Asymptotic behaviour of eigenfunctions for multipar¬ 
ticle schrodinger operators, Commun. Math. Phys., 34 (1973), pp. 251-270. 

[21] G. CsANYi, T. Albaret, G. Moras, M. Payne, and A. D. Vita, Multiscale hy¬ 
brid simulation methods for material systems, J. Phys,: Condens. Matter, 17 (2005), 
pp. 691-703. 

[22] G. CsANYi, T. Albaret, M. Payne, and A. D. Vita, “Learn on the fly”: a 
hybrid classical and quantum-mechanical molecular dynamics simulation, Phys. Rev. 
Lett., 93 (2004), pp. 175503 1-4. 

[23] W. E AND J. Lu, The elastic continuum limit of the tight binding model. Chin. Ann. 
Math. Ser. B, 28 (2007), pp. 665-675. 

[24] -, The electronic structure of smoothly deformed crystals: Cauchy-born rule for 

the nonlinear tight-binding model. Comm. Pure Appl. Math., 63 (2010), pp. 1432- 
1468. 

[25] -, The electronic structure of smoothly deformed crystals: Wannier functions and 

the cauchy-born rule. Arch. Ration. Mech. Anal., 199 (2011), pp. 407-433. 

[26] -, The Kohn-Sham equation for deformed crystals, Mem. Amer. Math. Soc., vol. 

221, no. 1040, 2013. 

[27] V. Ehrlacher, C. Ortner, and A. Shapeev, Analysis of boundary conditions 
for crystal defect atomistic simulations. arXiv: 1306.5334. 

[28] F. Ercolessi, Lecture notes on tight-binding molecular dynamics and tight-binding 
justification of classical potentials. Lecture notes 2005. 

[29] M. Finnis, Interatomic Forces in Condensed Matter, Oxford University Press, Ox¬ 
ford, 2003. 

[30] C. Fu AND K. Ho, First-principles calculation of the equilibrium ground-state prop¬ 
erties of transition metals: Applications to nb and mo, Phys. Rev. B, 28 (1983), 
pp. 5480-5486. 

[31] J. Gao and D. Truhlar, Quantum mechanical methods for enzyme kinetics, Annu. 
Rev. Phys. Chem., 53 (2002), pp. 467-505. 

[32] S. Goedecker, Integral representation of the fermi distribution and its applications 
in electronic-structure calculations, Phys. Rev. B, 48 (1993), pp. 17573-17575. 

[33] -, Linear scaling electronic structure methods. Rev. Mod. Phys., 71 (1999), 

pp. 1085-1123. 

[34] S. Goedecker and M. Teter, Tight-binding electronic-structure calculations and 
tight-binding molecular dynamics with localized orbitals, Phys. Rev. B, 51 (1995), 
pp. 9455-9464. 

[35] C. Goringe, D. Bowler, and E. Hernandez, Tight-binding modelling of ma¬ 
terials, Rep. Prog. Phys., 60 (1997), pp. 1447-1512. 


34 



[36] R. Horn and C. Johnson, Matrix Analysis, Cambridge University Press, Cam¬ 
bridge, 1991. 

[37] S. Ismail-Beigi and T. Arias, Locality of the density matrix in metals, semicon¬ 
ductors, and insulators, Phys. Rev. Lett., 82 (1999), pp. 2127-2130. 

[38] J. Kermode, T. Albaret, D. Sherman, N. Bernstein, P. Gumbsch, 
M. Payne, G. Gsanyi, and A. D. Vita, Low-speed fracture instabilities in a 
brittle crystal. Nature, 455 (2008), pp. 1224-1227. 

[39] W. Kohn, Analytic properties of block waves and wannier functions, Phys. Rev., 
115 (1959), pp. 809-821. 

[40] M. Kohyama and R. Yamamoto, Tight-binding study of grain boundaries in si: 
Energies and atomic structures of twist grain boundaries, Phys. Rev. B, 49 (1994), 
pp. 17102-17117. 

[41] G. Kresse and J. Furthmuller, Efficient iterative schemes for ab initio total- 
energy calculations using a plane-wave basis set, Phys. Rev. B, 54 (1996), pp. 11169- 
11186. 

[42] S. Lahbabi, The reduced Hartree-Eock model for short-range guantum crystals with 
nonlocal defects, Ann. Henri Poincare, 15 (2014), pp. 1403-1452. 

[43] S. Lee, J. Dow, and O. Sankey, Theory of charge-state splittings of deep levels, 
Phys. Rev. B, 31 (1985), pp. 3910-3914. 

[44] E. Lieb and B. Simon, The Thomas-Eermi theory of atoms, molecules and solids. 
Advances in Math., 23 (1977), pp. 22-116. 

[45] L. Lin, J. Lu, L. Ying, and W. E, Pole-based approximation of the fermi-dirac 
function. Chin. Ann. Math. B, 30 (2009), pp. 729-742. 

[46] M. Luskin and C. Ortner, Atomistic-to-continuum coupling, Acta Numerica, 
(2013). 

[47] P. Motamarria, M. Nowakb, K. Leiterg, J. Knapg, and V. Gavini, Higher- 
order adaptive finite-element methods for Kohn-Sham density functional theory, J. 
Comp. Phys., 253 (2013), pp. 308-343. 

[48] F. Nazar and C. Ortner. manuscript. 

[49] R. Nunes, J. Bennetto, and D. Vanderbilt, Structure, barriers, and relaxation 
mechanisms of kinks in the 90° partial dislocation in silicon, Phys. Rev. Lett., 77 
(1996), pp. 1516-1519. 

[50] S. Ogata, E. Lidorikis, F. Shimojo, A. Nakano, P. Vashishta, and 
R. Kalia, Hybrid finite-element/molecular-dynamic/electronic-density-functional 
approach to materials simulations on parallel computers, Comput. Phys. Commun., 
138 (2001), pp. 143-154. 


35 



[51] C. Ortner and F. Theil, Justification of the cauchy-born approximation of elas- 
todynamics, Arch. Ration. Mech. Anal., 207 (2013), pp. 1025-1073. 

[52] D. Papaconstantopoulos, Handbook of the Band Structure of Elemental Solids, 
From Z = 1 To Z = 112, Springer New York, 2015. 

[53] E. Prodan and W. Kohn, Nearsightedness of electronic matter, Proc. Natl. Acad. 
Sci. USA, 102 (2005), pp. 11635-11638. 

[54] J. Slater and G. Koster, Simplified Icao method for the periodic potential prob¬ 
lem, Phys. Rev., 94 (1954), pp. 1498-1524. 

[55] J. Tersoff, New empirical approach for the structure and energy of covalent sys¬ 
tems, Phys. Rev. B, 37 (1988), pp. 6991-7000. 

[56] C. Wang, C. Chan, and K. Ho, Tight-binding molecular-dynamics study of defects 
in silicon, Phys. Rev. Lett., 66 (1991), pp. 189-192. 

[57] Y. Wang, G. Stogks, W. Shelton, D. Nigholson, Z. Szotek, and W. Tem- 
MERMAN, Order-N multiple scattering approach to electronic structure calculations, 
Phys. Rev. Lett., 75 (1995), pp. 2867-2870. 

[58] A. Warshela and M. Levitta, Theoretical studies of enzymic reactions: Di¬ 
electric, electrostatic and steric stabilization of the carbonium ion in the reaction of 
lysozyme. Journal of Molecular Biology, 103 (1976), pp. 227-249. 

[59] H. Ylfkawa, On the interaction of elementary particles, Proc. Phys. Math. Soc. 
Japan., 17 (1935), pp. 48-57. 


36 



